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Abstract 

Fluctuating liydrodynamics is used to describe the total energy fluctuations of a freely evolving 
gas of inelastic hard spheres near the threshold of the clustering instability. They are shown 
to be governed by vorticity fluctuations only, that also lead to a renormalization of the average 
total energy. The theory predicts a power-law divergent behavior of the scaled second moment of 
the fluctuations, and a scaling property of their probability distribution, both in agreement with 
simulations results. A more quantitative comparison between theory and simulation for the critical 
amplitudes and the form of the scaling function is also carried out. 
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Granular gases are assemblies of macroscopic particles evolving independently between 
inelastic collisions Q]. The methods of non-equilibrium statistical mechanics, kinetic theory, 
and hydrodynamics have been successfully extended to describe the observed macroscopic 
behavior and also, although in a much more limited form, the fluctuations around it O]. The 
lack of energy conservation makes them to behave quite differently from molecular fluids. A 
simple widely used model for granular gases consists of smooth inelastic hard spheres (IHS), 
with momentum conserving dynamics. Inelasticity is characterized by means of a constant 
coefficient of normal restitution a. 

Recently, molecular dynamics (MD) simulation results have been reported for the total 
energy fluctuations of a two-dimensional freely evolving IHS gas, near the threshold of the 
clustering instability The dimensionless second moment was found to exhibit a power- 
law divergent behavior with the distance to the instability. Also, the scaled cooling rate 
was found to tend to zero according to a power law, although in a weak way. Besides, 
the distribution function for the energy fluctuations, when properly scaled, turned out to 
be independent of the parameters defining the system. This was associated with a scaling 
property of the distribution. Quite remarkably, the scaling function was very well fitted by 



the same expression as several equilibrium and non-equilibrium molecular systems 



The main goal of the present Letter is to provide an explanation for the above results on 

n 

the basis of fluctuating hydrodynamics ;6], showing that it gives an accurate description of 
the fluctuations occurring in a granular fluid in the threshold of the clustering instability. 

Consider an isolated system of IHS of mass m and diameter o. The total (kinetic) 
energy E of the system can be expressed in the form 

E{t) = ^Jdr [dn{r, t)f{r, t) + mn{r, t)u^{r, t)] , (1) 

where d is the dimension of the system, n{r,t) the number density field, T{r,t) the tem- 
perature field, and u{r,t) the flow field. The tildes indicate that all the quantities are 
understood as fluctuating variables. The above expression can be justified by identifying 
the definitions of the microscopic densities with their fluctuating values. In the following, 
systems in the homogeneous cooling state (HCS) will be considered. At a macroscopic level, 
this state is characterized by a constant uniform density uh, a vanishing flow field uh = 0, 
and a uniform time-dependent temperature obeying the law (if dtTnit) = —CH(TH)TH(t), 
where TnitY^"^ is the cooling rate. Moreover, we will restrict ourselves to the region in 
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which the amphtudes of the fluctuations of the flelds around their HCS values remain small 
on the average. Then retaining up to quadratic order in the deviations, Eq. (Q) yields 



5Eit) = E{t)-EH{t) 

dnnSfir, t) + d6h{r, t)6f{r, t) 



- / dr 
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+ mnH\5u{r,t)\'^]. (2) 

Here, Euif) = dNTH{t)/2, 6h{r,t) = n{r,t)-nH, 6u{r,t) = u{r,t), and 5T(r,t) = f{r,t)- 
Tnit). It is now convenient to introduce dimensionless position, I, and time, s, scales as 
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I = v/Iq and ds = VHit)dt/lQ, respectively, where vh = [2T//(t)/m] ' is the termal velocity 
and Iq = {nH(j'^~^)~^ is proportional to the mean free path. Moreover, dimensionless flelds 
are deflned by p{l,s) = 6h{r,t)/nH, <^{l,s) = 5u{r,t)/vH{t), and 9{l,s) = 5T{r,t)/TH{t). 
Then, Eq. Q takes the form 



k 
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Pk{s)9^k{s) + ^l^fc(s)| 



(3) 



where e(s) = 6E{s)/Eh{s), V = L'^ is the volume of the system in the new units, and 
the Fourier transforms of the flelds have been introduced. It is assumed that after a time 
of the order of the mean free time between collisions, the system reaches a hydrodynamic 
regime in which all the energy of the system is stored in the hydrodynamic modes. In 
this regime, the hydrodynamic flelds are expected to be described at a mesoscopic level 
by fluctuating hydrodynamic equations. Here, they will be assumed to be linear Langevin 
equations obtained by linearizing the Navier-Stokes equations for a granular gas around 
the HCS. Moreover, the assumption is made that the noise terms are deflned by the same 
properties as for molecular, elastic gases P]. This is not expected to be true, except in 
the nearly elastic limit, i.e. when the coefficient of normal restitution a is very close to 
unity. Consequently, the theory will be restricted in the following to this limit. Thus, the 



transversal flow fleld or vorticity fleld, uJk±, obeys the equation 



{ds - C72 + V*k') ujk±is) = (4) 

In this expression, (* = CH[TH(t)]io/vH{t) and rj* = riH[TH{t)]/mnHloVH(t), rjn being the 
shear viscosity. The random noise term ^k±(s) is Gaussian, with a correlation 

{^k±{s)^k'±{s')) = ^6(3 - s')5k,-k'V*k\ (5) 



I being the unit tensor of dimension d — 1 in the subspace perpendicular to k, and the 
angular brackets denoting average over the noise realizations. A main advantage of using 
the scaled variables is that the coefficients in Eq. (0]) and the strength of the noise become 
time-independent, contrary to what happens in the original variables. The equation shows 
that uJkj_ grows in time for those values of k such that X±{k) = C*/2 — 1]*^^"^ > 0- Although 
this does not imply by itself that the HCS is linearly unstable, due to the time- dependent 
scaling of the velocity introduced above, simulation results and nonlinear analytical analysis 
of the Navier-Stokes equations have shown that this growth is the origin of the clustering 
instability The minimum value of k for a system of linear extent L, measured in the 

/-scale, is kmm = 2tt/L. Then, for given values of the other parameters, the system becomes 
unstable if L > L^, with = 27r(2?7*/C*)^/^. For L < L^, the HCS is stable and the long 
time solution of Eq. (jH) is 

u^kL{s)= r ds'e^^-^'^'^^'^^^^is'). (6) 
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From this expression, it is easily obtained 

for s > s' ^ 1. This shows that as L tends to Lc from below, the amplitudes of the 
fluctuations of the transversal components of the velocity increase very fast as the instability 
is approached due to contributions from values of k close to kc- For the same reason the 
decay of these fluctuations becomes very slow. This is not the case for the fluctuations 
of the other hydrodynamic fields, whose Langevin equations are decoupled from Eq. 

Therefore, it seems possible to consider a range of values of 6L = {Lc — L)/Lc where 
the fiuctuations of uJk,±, although still small, dominate over the fiuctuations of density and 
temperature. But, although this is true for components with A; > 0, some care is needed 
when analyzing Eq. Q, since it involves the A; = component of the temperature field, 
9o{s). The Langevin equation for e(s) is obtained from the linearization around the HCS of 
the macroscopic average equation for the total energy, 

dtE{t) = J drn{r,t)CH{n,T)T{r,t). (8) 

The result is: 



(9) 



Here, the dependence of the coohng rate on the temperature has been taken into account. 



Moreover, the noise term discussed in [1^, associated with the locahzed character of the 
energy dissipation, has been omitted. Ahhough it can be expected to be neghgible far 
from the instabihty in the quasi-elastic hmit, this may be not the case near the instabihty. 
Equation Q shows a couphng between the fluctuations of the volume averaged temperature 
and those of the total energy. Use of Eq. Q into Eq. Q and neglecting contributions from 



the density and longitudinal velocity fluctuations gives 



UJ[S 



J2 1'^fc^' 



(10) 



valid in the region 6L <^1. The long time limit of the average value of e(s) is, therefore 

3{d 



{e)st = lim {uj{s)) 



Nd 



1) ^ r]*e 



(11) 



Since we are considering -C 1, the sum over k in the above expression is dominated by 
the 2d modes with the largest wavelength, for which X±{kmin) — ~C*SL. Using this into Eq. 
(fTTj) . it follows that there is a renormalization by fluctuations of the average total energy of 
the HCS, E{t) = {E{t)), given by 

3{d-l) 



E{t) = Enit) 



1 + 



6L 



(12) 



Similarly, there is also a renormalization of the temperature of the HCS, T{t) = {T[t))st, 
that can be evaluated directly from the long time limit of the average of Eq. (jUj), 



Tit) = Tnit) 



V 



Tnit) 



1 + 



(13) 



Alternatively, an effective temperature T^fit) can be defined as T^fif) = 2E(t)/Nd. Of 
course, the form of the renormalized law for the temperature depends on the definition used 
for the latter. In what was actually measured was Ce/ = Cef(Tef)lo/vH(Tef), with (ef 
defined by dtTej = —(ef(Tef)Tef. Then, it is found 



3(rf-l)~-i 



(14) 



This result predicts that near the clustering instability threshold, (C, 
Af^6L with = 3{d — l)/n//L^, that is just the behavior observed in 



*-2 



*-2 
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Define now SE{s) = [E{s) — E{s)]/E{s) = [e{s) — {e)st]EHis)/E{s). We are considering 
deviations from tlie renormalized average value, i.e. including the fluctuations effects, and 
not from the macroscopic bare value. A standard calculation using Eq. (0) and exploiting 
the Gaussian character of the noise, gives that in the instability threshold and for s > s' ^ 1 
it is 

(Pis) - m^ms') - mst])st = -fj^ , (15) 

where Sc = {2(*6L)~^ is a divergent "critical" relaxation time. Now Eq. can be easily 
solved with the result, when 6L <^ 1, 

{6E{s)6E{s'))st = {p{s) - {oJ)stMs') - {iJ)st])st, (16) 

valid for s > s' ^ 1. Thus below the instability, the scaled total energy fluctuations decay 
with the same rate as the fluctuations of the kinetic energy associated with the transversal 
modes of the velocity. For s = s', Eq. ()16p yields 

4 ^ {m')st = AirL\ (17) 

with = 9((i— Vj/n'^jL'jf'd. Therefore, close to the instability point, the relative dispersion 
of the total energy fluctuations (Te presents a divergent behavior with a critical exponent —1, 
and an amplitude depending on the density uh and the coefficient of normal restitution 
a (through the value of the critical length L^). Again, this is the same behavior as reported 
in 0] from MD simulations. 

To carry out a more detailed check of the theory presented here, we have performed MD 
simulations of two-dimensional systems with different values of a and Tin (see Table In all 
cases, the dependence on 5L of both the cooling rate and the dispersion of the total energy, 
i.e. the exponents in the power laws (jl4j) and (jl7|) . was in agreement with the theoretical 
predictions. This was illustrated in Figs. 1 and 2 of ref. P| and no more details will be given 
here. The comparison between the predicted critical amplitudes and the MD results given 
in Table |l] can be considered as satisfactory, in the sense that the theory correctly predicts 
the order of magnitude of the amplitudes, specially taking into account the smallness of the 
quantities being measured. 

Next, let us proceed to investigate the form of the probability distribution of the energy 
fluctuations. Particularization of Eq. for the modes with the smallest possible value of 
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TABLE I: Comparison between the predicted and MD values for the critical amplitudes of the 
cooling rate, A(^, and the total energy dispersion, A^. All the values of the amplitudes have been 
multiplied by 10'^. 



0.02 0.9 0.88 1.11 0.62 0.6 

0.02 0.8 1.62 3.63 1.15 1.5 

0.1 0.98 1.06 0.50 0.75 0.47 

0.1 0.95 2.4 2.38 1.7 1.45 

0.2 0.98 1.97 2.34 1.4 1.34 

0.2 0.95 4.59 7.78 3.24 3.6 



k in the limit 6L <^ 1 gives 

{ds + CSL)u;k±{s) = ^k^{s), (18) 

where it is understood that \k\ = kmin- Define a new time scale dr = (*SL ds, and a new 
transversal velocity field by cl'^_|_ = cjfcx/Lfcr]/^. Equation (fTHj) becomes 

{d^ + l)ul^ = a^{T), (19) 

with 

{aAr)a'Ar')) = ^'^^-"'^(^ " (20) 

Equations (fT^ and (pUj) imply that the probability distribution for uj1.j_, with |fc| = fcmm, 
near the clustering instability depends only on the dimension d of the system. In fact, since 
the noise term ^^^(t) is Gaussian, it is trivial to write the long time form of this distribution 
using Eq. ((Tj) with s = s', 

p,.Kj = (2..y-'^-"'%-"a/-J, (21) 

with 0"^ = (i^/^/12((i— 1)^/^. In the time scale r, and keeping only the dominant modes, Eq. 
(fTn|) reads 

SLdry = -l(y--^ E \^lAr)n, (22) 



2 r d . , , 



mm 



where y = e/as and the sum is restricted to vectors k with |fc| = fcmm- From the comparison 
of Eqs. (fT^ and it is seen that, on the r scale and in the threshold of the instability, 



y decays much faster than the dominant components of w^j^. Consequently, for large r 
the solution of Eq. is y = §Z]|fc|=fc„i„ l'*^fc±(''")P) where the probability distributions of 
the modes uj^^^ is given by Eq. Since the later does not depend on the parameters 

of the system other than the dimensionality, the same property follows for the probability 
distribution of both y and the variable 

^-l = -[d{d-i)Y'^ + l E KJ'- (23) 



(Te d 



k\=k 



mi 71 



This is equivalent to saying that the probability distribution for 6E verifies the scaling 
relation 



pm = -/ ( ^ 1 , (24) 



f5E' 
cte \<^e, 

where / is a scaling function. This is just the property assumed in j2j and verified by MD 
simulations. Since the probability distribution function for cj^j^ is known, it is possible to 
numerically generate the probability distribution function P{6E). The result is shown in 
Fig. HJ Also plotted is the function 

n(5E) = ir(e^-^")^ x = -b{s + 6E), a = 7r/2, (25) 

with K = 2.U,b = 0.938, and s = 0.374, that fits extremely well the MD results for aEP(5E) 
P]. It is important to remark that fluctuations in a large number of equilibrium and non- 



equilibrium systems exhibiting self-organized criticality as well as confined turbulent flows, 
present the same kind of behavior [J, 0]. Although the agreement between both plotted 
curves is not so bad for positive values of 5E, strong discrepancies are observed for negative 
values. A major source for them is easily identified from Eq. (P^ . that for d = 2 implies 
5E/(Je > — "\/2, while smaller values are found in the MD simulations. Since Eq. fj23p is a 
consequence of Eq. 0, it seems plausible that in order to elaborate a more accurate theory 
the intrinsic noise associated with the cooling rate must be taken into account. 

In summary, we have developed a mesoscopic theory for the fluctuations of the total 
energy of an isolated granular gas near the threshold of the clustering instability. The the- 
ory describes accurately the qualitative behavior obtained in MD simulations, namely the 
divergent behavior of the dimensionless second moment and the decrease of the apparent 
cooling rate. Also, it is consistent with the observed scaling property of the probability dis- 
tribution function of the fluctuations. On the other hand, it seems clear that a more refined 
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FIG. 1: Probability density function of the relative total energy fluctuations aEPhi^E) for a system 
of inelastic hard disks. The broken line is the theoretical prediction derived in this paper and the 
solid line is Eq. 

formulation is needed in order to get a more satisfactory quantitative agreement, especially 
for the distribution function. In any case, we believe the present work clearly indicates the 
way in which fluctuations in granular systems near an instability can be analyzed and, in 
particular, the dominant role played by nonlinear coupling between hydrodynamic modes. 
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